Modele liniowe

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns

Na początek ustawimy globalnie opcję set_output ze Scikit-learn'a (dokumentacja). Dzięki temu transformowane dane będą zawsze na wyjściu typu pd.DataFrame, zamiast np.ndarray, i zachowają odpowiednie nazwy cech.

import sklearn

sklearn.set_config(transform_output="pandas")

Przewidywanie cen domów

Wykorzystamy zbiór danych Ames housing, w którym zadaniem jest przewidywanie wartości domu na podstawie cech budynku, działki, lokalizacji itp. Jest to więc przewidywanie wartości ciągłej, czyli regresja. Obejmuje zmienne do usunięcia (np. ID transakcji), numeryczne (floaty i inty), kategoryczne nieuporządkowane (categorical nominal) oraz kategoryczne uporządkowane (categorical ordinal), więc będzie wymagał wstępnego przetworzenia tak jak większość prawdziwych danych w ML.

Inne znane, ale gorsze jakościowo zbiory tego typu to na przykład:

Autor zbioru to Dean De Cock, a zbiór został opisany oryginalnie w tym artykule. Szczegółowe opisy zmiennych znajdują się w pliku ames_description.txt.

Formatem pliku jest Apache Parquet, bardzo wydajny format binarny. Ma on szereg zalet w ML i inżynierii danych porównaniu do plików CSV:

Dla zainteresowanych szczegółami: link 1, link 2, link 3. Ten format ma szerokie wsparcie, w tym w Pandasie. W przypadku tego zbioru, CSV zajmuje ~1.2 MB, a Parquet ~0.2 MB.

df = pd.read_parquet("ames_data.parquet")

# remove dots from names to match data_description.txt
df.columns = [col.replace(".", "") for col in df.columns]

df.head()
Order PID MSSubClass MSZoning LotFrontage LotArea Street Alley LotShape LandContour ... PoolArea PoolQC Fence MiscFeature MiscVal MoSold YrSold SaleType SaleCondition SalePrice
0 1 526301100 20 RL 141.0 31770 Pave None IR1 Lvl ... 0 None None None 0 5 2010 WD Normal 215000
1 2 526350040 20 RH 80.0 11622 Pave None Reg Lvl ... 0 None MnPrv None 0 6 2010 WD Normal 105000
2 3 526351010 20 RL 81.0 14267 Pave None IR1 Lvl ... 0 None None Gar2 12500 6 2010 WD Normal 172000
3 4 526353030 20 RL 93.0 11160 Pave None Reg Lvl ... 0 None None None 0 4 2010 WD Normal 244000
4 5 527105010 60 RL 74.0 13830 Pave None IR1 Lvl ... 0 None MnPrv None 0 3 2010 WD Normal 189900

5 rows × 82 columns

df.shape
(2930, 82)

Eksploracja danych, czyszczenie danych i inżynieria cech

Wstępne czyszczenie danych (data cleaning) obejmuje:

To drugie jest motywowane wykresem przedstawionym poniżej. Zostało uznane za błąd już przez samego autora.

df = df.drop(["Order", "PID"], axis="columns")
df = df.loc[~df["Neighborhood"].isin(["GrnHill", "Landmrk"]), :]

plt.scatter(df["GrLivArea"], df["SalePrice"])
plt.title("House area vs price")
plt.xlabel("GrLivArea")
plt.ylabel("SalePrice")
plt.show()

df = df.loc[df["GrLivArea"] <= 4000, :]

plt.scatter(df["GrLivArea"], df["SalePrice"])
plt.title("House area vs price, outliers removed")
plt.xlabel("GrLivArea")
plt.ylabel("SalePrice")
plt.show()

Zawsze warto też przyjrzeć się rozkładowi zmiennej docelowej, żeby poznać jego typ i skalę. Jak widać poniżej, rozkład jest dość skośny, co ma sens - mało jest bardzo drogich domów. Można to zatem uznać za problem regresji niezbalansowanej (imbalanced regression), w którym rozkład zmiennej docelowej nie jest normalny.

Załóżmy, że typowego klienta nie interesują bardzo drogie domy. Cały ciężki ogon rozkładu można zatem uznać zasadniczo za outliery - im większa wartość, tym mniej interesuje przeciętną osobę. Warto też pamiętać, że realny "przeciętny" dochód to mediana, a nie średnia zarobków, co odpowiednio przekłada się na zainteresowanie domami.

W tym przypadku warto przetransformować tę zmienną logarytmem bliżej rozkładu normalnego z dwóch powodów.

Po pierwsze, stabilność numeryczna - wartości są bardzo duże, rzędu setek tysięcy, i przy obliczaniu odchylenia standardowego czy wariancji to byłby problem. Logarytm to przykład transformacji stabilizującej wariancję (variance-stabilizing transform).

Ponadto regresja liniowa (pierwszy model liniowy, który będziemy rozważać) czyni szereg założeń w modelowaniu, co wynika z jej podstaw statystycznych. Są to między innymi:

  1. Normalność błędów (residuals normality) - błędy modelu, czyli y − , mają rozkład normalny
  2. Homoskedastyczność (homoscedasticity) - błędy mają stałą wariancję, czyli nie zależą od wielkości zmiennej zależnej

W naszym przypadku oznacza to, że podobnie często nie doceniamy albo przeceniamy wartość domu, oraz dla tanich i drogich domów mamy podobny błąd. Przy rozkładzie zmiennej zależnej bardzo różnym od normalnego te założenia są typowo łamane, co skutkuje gorszą jakością predykcji modelu. Dlatego transformacji warto też dokonać dla lepszej jakości modelu.

Skąd wiemy, że logarytm? Bo rozkłady o dodaniej skośności (positively skewed) przesuwa się "w prawo" właśnie logarytmem, a w drugą stronę eksponentą. Dla zainteresowanych: link 1, link 2. Są też algorytmy, które potrafią same estymować potrzebną transformację, np. Box-Cox transform albo Yeo-Johnson transform. Podobne transformacje można też oczywiście wykonywać nie tylko dla zmiennej zależnej, ale też (a nawet jest to częstsze) dla cech.

Zadanie 1 (1 punkt)

  1. Narysuj wykres cen - histogram kolumny "SalePrice". Pamiętaj o podaniu tytułu i dodaniu legendy.
  2. Oblicz średnią i medianę cen, a następnie nanieś je na wykres jako pionowe, przerywane linie. Może się przydać ten przykład. Wypisz te wartości, lub zawrzyj je w opisie w legendzie, zaokrąglone do liczb całkowitych.
  3. Zlogarytmuj ceny z pomocą np.log1p (zamień istniejące wartości).
  4. Narysuj wykres po zlogarytmowaniu, razem z liniami dla średniej i mediany.

Może ci się przydać metoda .plot.hist() dla pd.Series. Zwraca ona normalny wykres Matplotliba, do którego można potem dodawać kolejne elementy.

# Plot before taking logarithm
plt.figure(figsize=(10, 6))
df["SalePrice"].plot.hist(edgecolor='black', bins=30, alpha=0.7, label="Sale Price")
plt.title("Sale Price Histogram")
plt.xlabel("Sale Price")
plt.ylabel("Number of holders")
plt.legend()

mean_price = np.mean(df["SalePrice"])
median_price = np.median(df["SalePrice"])

print(f"Średnia: {int(mean_price)}")
print(f"Mediana: {int(median_price)}")

plt.axvline(mean_price, color='red', linestyle="--")
plt.axvline(median_price, color='blue', linestyle="--")
plt.text(mean_price, 0 , "mean", color="red", rotation=60, fontsize = 12)
plt.text(median_price, 100 , "median", color="blue", rotation=60, fontsize = 12)

# Plot after taking logarithm
df["SalePrice"] = np.log1p(df["SalePrice"])
plt.figure(figsize=(10, 6))
df["SalePrice"].plot.hist(edgecolor='black', bins=30, alpha=0.7, label="Sale Price")
plt.title("Sale Price Histogram after taking logarithm")
plt.xlabel("Sale Price")
plt.ylabel("Number of holders")
plt.legend()

mean_price_log = np.mean(df["SalePrice"])
median_price_log = np.median(df["SalePrice"])

print(f"Średnia (po zlogarytmowaniu): {int(mean_price_log)}")
print(f"Mediana (po zlogarytmowaniu): {int(median_price_log)}")

plt.axvline(mean_price_log, color='red', linestyle='--')
plt.axvline(median_price_log, color='blue', linestyle='--')
plt.text(mean_price_log, 0 , "mean", color="red", rotation=60, fontsize = 12)
plt.text(median_price_log, 100 , "median", color="blue", rotation=60, fontsize = 12)

plt.show()
Średnia: 180358
Mediana: 160000
Średnia (po zlogarytmowaniu): 12
Mediana (po zlogarytmowaniu): 11

W zbiorze znajdują się liczne wartości brakujące, które zostały uzupełnione na podstawie tego notebooka na Kaggle.

def replace_na(df: pd.DataFrame, col: str, value) -> None:
    df.loc[:, col] = df.loc[:, col].fillna(value)

# Alley : data description says NA means "no alley access"
replace_na(df, "Alley", value="None")

# BedroomAbvGr : NA most likely means 0
replace_na(df, "BedroomAbvGr", value=0)

# BsmtQual etc : data description says NA for basement features is "no basement"
replace_na(df, "BsmtQual", value="No")
replace_na(df, "BsmtCond", value="No")
replace_na(df, "BsmtExposure", value="No")
replace_na(df, "BsmtFinType1", value="No")
replace_na(df, "BsmtFinType2", value="No")
replace_na(df, "BsmtFullBath", value=0)
replace_na(df, "BsmtHalfBath", value=0)
replace_na(df, "BsmtUnfSF", value=0)

# Condition : NA most likely means Normal
replace_na(df, "Condition1", value="Norm")
replace_na(df, "Condition2", value="Norm")

# External stuff : NA most likely means average
replace_na(df, "ExterCond", value="TA")
replace_na(df, "ExterQual", value="TA")

# Fence : data description says NA means "no fence"
replace_na(df, "Fence", value="No")

# Functional : data description says NA means typical
replace_na(df, "Functional", value="Typ")

# GarageType etc : data description says NA for garage features is "no garage"
replace_na(df, "GarageType", value="No")
replace_na(df, "GarageFinish", value="No")
replace_na(df, "GarageQual", value="No")
replace_na(df, "GarageCond", value="No")
replace_na(df, "GarageArea", value=0)
replace_na(df, "GarageCars", value=0)

# HalfBath : NA most likely means no half baths above grade
replace_na(df, "HalfBath", value=0)

# HeatingQC : NA most likely means typical
replace_na(df, "HeatingQC", value="Ta")

# KitchenAbvGr : NA most likely means 0
replace_na(df, "KitchenAbvGr", value=0)

# KitchenQual : NA most likely means typical
replace_na(df, "KitchenQual", value="TA")

# LotFrontage : NA most likely means no lot frontage
replace_na(df, "LotFrontage", value=0)

# LotShape : NA most likely means regular
replace_na(df, "LotShape", value="Reg")

# MasVnrType : NA most likely means no veneer
replace_na(df, "MasVnrType", value="None")
replace_na(df, "MasVnrArea", value=0)

# MiscFeature : data description says NA means "no misc feature"
replace_na(df, "MiscFeature", value="No")
replace_na(df, "MiscVal", value=0)

# OpenPorchSF : NA most likely means no open porch
replace_na(df, "OpenPorchSF", value=0)

# PavedDrive : NA most likely means not paved
replace_na(df, "PavedDrive", value="N")

# PoolQC : data description says NA means "no pool"
replace_na(df, "PoolQC", value="No")
replace_na(df, "PoolArea", value=0)

# SaleCondition : NA most likely means normal sale
replace_na(df, "SaleCondition", value="Normal")

# ScreenPorch : NA most likely means no screen porch
replace_na(df, "ScreenPorch", value=0)

# TotRmsAbvGrd : NA most likely means 0
replace_na(df, "TotRmsAbvGrd", value=0)

# Utilities : NA most likely means all public utilities
replace_na(df, "Utilities", value="AllPub")

# WoodDeckSF : NA most likely means no wood deck
replace_na(df, "WoodDeckSF", value=0)

# CentralAir : NA most likely means No
replace_na(df, "CentralAir", value="N")

# EnclosedPorch : NA most likely means no enclosed porch
replace_na(df, "EnclosedPorch", value=0)

# FireplaceQu : data description says NA means "no fireplace"
replace_na(df, "FireplaceQu", value="No")
replace_na(df, "Fireplaces", value=0)

# SaleCondition : NA most likely means normal sale
replace_na(df, "SaleCondition", value="Normal")

# Electrical : NA most likely means standard circuit & breakers
replace_na(df, "Electrical", value="SBrkr")

W dalszej kolejności zamienimy zmienne MSSubClass i MoSold z numerycznych na kategoryczne, zgodnie z ich znaczeniem. Zakodujemy także zmienne kategoryczne uporządkowane (categorical ordinal) z tekstowych na kolejne liczby całkowite.

df = df.replace(
    {
        "MSSubClass": {
            20: "SC20",
            30: "SC30",
            40: "SC40",
            45: "SC45",
            50: "SC50",
            60: "SC60",
            70: "SC70",
            75: "SC75",
            80: "SC80",
            85: "SC85",
            90: "SC90",
            120: "SC120",
            150: "SC150",
            160: "SC160",
            180: "SC180",
            190: "SC190",
        },
        "MoSold": {
            1: "Jan",
            2: "Feb",
            3: "Mar",
            4: "Apr",
            5: "May",
            6: "Jun",
            7: "Jul",
            8: "Aug",
            9: "Sep",
            10: "Oct",
            11: "Nov",
            12: "Dec",
        },
    }
)
df = df.replace(
    {
        "Alley": {"None": 0, "Grvl": 1, "Pave": 2},
        "BsmtCond": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "BsmtExposure": {"No": 0, "Mn": 1, "Av": 2, "Gd": 3},
        "BsmtFinType1": {
            "No": 0,
            "Unf": 1,
            "LwQ": 2,
            "Rec": 3,
            "BLQ": 4,
            "ALQ": 5,
            "GLQ": 6,
        },
        "BsmtFinType2": {
            "No": 0,
            "Unf": 1,
            "LwQ": 2,
            "Rec": 3,
            "BLQ": 4,
            "ALQ": 5,
            "GLQ": 6,
        },
        "BsmtQual": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "ExterCond": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "ExterQual": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "FireplaceQu": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "Functional": {
            "Sal": 1,
            "Sev": 2,
            "Maj2": 3,
            "Maj1": 4,
            "Mod": 5,
            "Min2": 6,
            "Min1": 7,
            "Typ": 8,
        },
        "GarageCond": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "GarageQual": {"No": 0, "Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "HeatingQC": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "KitchenQual": {"Po": 1, "Fa": 2, "TA": 3, "Gd": 4, "Ex": 5},
        "LandSlope": {"Sev": 1, "Mod": 2, "Gtl": 3},
        "LotShape": {"IR3": 1, "IR2": 2, "IR1": 3, "Reg": 4},
        "PavedDrive": {"N": 0, "P": 1, "Y": 2},
        "PoolQC": {"No": 0, "Fa": 1, "TA": 2, "Gd": 3, "Ex": 4},
        "Street": {"Grvl": 1, "Pave": 2},
        "Utilities": {"ELO": 1, "NoSeWa": 2, "NoSewr": 3, "AllPub": 4},
    }
)

W Pandasie zmienne kategoryczne, napisy itp. są typu object, natomiast numeryczne mają typy z Numpy'a, np. int64. Oczywiście oba rodzaje zmiennych trzeba przetwarzać na inne sposoby, więc zapiszemy to od razu. Wyodrębnimy też wektor docelowy.

Od razu dokonamy też podziału na zbiór treningowy i testowy. Jako że nasz zbiór jest dość mały, to podział będzie 70%-30%.

from sklearn.model_selection import train_test_split

y = df.pop("SalePrice")

categorical_features = df.select_dtypes(include="object").columns
numerical_features = df.select_dtypes(exclude="object").columns

X_train_raw, X_test_raw, y_train, y_test = train_test_split(
    df, y, test_size=0.3, random_state=0
)

Analiza i przetwarzanie cech

Modele liniowe działają najlepiej, kiedy nasze cechy są odpowiednio przetworzone:

Na start imputujemy wartości brakujące w zmiennych numerycznych i przeskalujemy je do zakresu [0,1]. Potem dokonamy analizy, czy mamy mocno skorelowane zmienne.

Przy przetwarzaniu danych krok po kroku lepiej robić nowe zmienne, żeby dało się odpalić znowu komórkę notebooka. Przy reużywaniu zmiennych mogłaby nastąpić konieczność uruchomienia wszystkiego od początku. Może ci się to przydać w zadaniach.

Pamiętaj, żeby we wszystkich kolejnych zadaniach zastosować taką samą transformację na zbiorze treningowym i testowym. Wszelkie wartości estymujemy na zbiorze treningowym.

Zadanie 2 (0.5 punktu)

Wybierz same zmienne numeryczne. Dokonaj imputacji wartości brakujących medianą, oraz przeskaluj wartości do zakresu [0,1]. Przydatne klasy: ColumnTransformer, Pipeline, SimpleImputer, MinMaxScaler.

Domyślnie ColumnTransformer dopisuje do nazw cech wyjściowych nazwę transformera jako prefix. Bywa to dość nieczytelne i problematyczne w dalszym przetwarzaniu danych. Żeby to wyłączyć, przekaż opcję verbose_feature_names_out=False.

from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import MinMaxScaler

median_imputer = SimpleImputer(strategy='median')
min_max_scaler = MinMaxScaler()

numerical_transformer = Pipeline(steps=[
    ('imputer', median_imputer),
    ('scaler', min_max_scaler)
])

preprocessor = ColumnTransformer(
    transformers=[
        ('num', numerical_transformer, numerical_features)
    ], verbose_feature_names_out=False
)

X_train_trans = preprocessor.fit_transform(X_train_raw)
X_test_trans = preprocessor.transform(X_test_raw)

Sprawdźmy teraz, czy mamy mocno skorelowane zmienne. Typowo przyjmuje się wartość bezwzględną korelacji powyżej progu 0.8/0.9/0.95. Pierwszą weryfikację można zrobić "na oko" na podstawie wykresu typu heatmap. Przekątną oczywiście ignorujemy - każda cecha ma korelację 1.0 z samą sobą.

Jeżeli zauważymy grupy skorelowanych cech, to z każdej grupy trzeba zostawić jedną cechę. Biblioteka Feature-engine, komplatybilna ze Scikie-learn, oferuje tutaj wygodną klasę SmartCorrelatedSelection, która pozwala z każdej grupy cech zostawić tę o największej wariancji.

Czemu o największej wariancji? Bo taka cecha jest mocno zmienna, a więc niesie jakąś informację. Cechy o bardzo niskiej wariancji praktycznie się nie zmieniają, więc prawdopodobnie nie są zbyt informatywne.

Zadanie 3 (0.5 punktu)

Oblicz korelacje między cechami w zbiorze treningowym, oraz narysuj heatmapę korelacji. Pandas potrafi obliczyć korelacje, a Seaborn - narysować wygodnie heatmapę.

Następnie wykorzystaj Feature-engine do eliminacji mocno skorelowanych cech, z progiem 0.8 (mamy sporo cech, więc eliminujmy, czemu nie), pozostawiając cechy o największej wariancji. Wypisz nazwy wyeliminowanych cech.

correlation_matrix = X_train_trans.corr()

plt.figure(figsize=(20, 16))
sns.heatmap(correlation_matrix, cmap="turbo")
plt.title("Features correlations heatmap")
plt.show()

from feature_engine.selection import SmartCorrelatedSelection

smart_correlated_selector = SmartCorrelatedSelection(
    threshold=0.8,
    selection_method="variance"
)

smart_correlated_selector.fit(X_train_trans, y_train)
X_train_num = smart_correlated_selector.transform(X_train_trans)
X_test_num = smart_correlated_selector.transform(X_test_trans)

eliminated_features = smart_correlated_selector.features_to_drop_
print("Eliminated features: ", eliminated_features)
Eliminated features:  ['BsmtFinSF2', 'TotRmsAbvGrd', 'Fireplaces', 'GarageArea', 'GarageCond', 'PoolQC']

Przejdźmy teraz do zmiennych kategorycznych. Jeżeli pewne kategorie występują bardzo rzadko, to prawdopodobnie nie ma sensu ich rozróżniać, i można by je połączyć w jedną, np. "Other" albo "Rare". Eliminacja rzadkich kategorii bardzo zmniejszy wymiarowość po one-hot encodingu, redukując złożoność obliczeniową i pamięciową. Często robi się to nawet głównie z tego powodu, nawet jeżeli lekko pogorszy to wyniki (aczkolwiek często daje nieco lepsze).

Po tym kroku możemy już wykonać one-hot encoding i uzyskać nasze zmienne kategoryczne.

Zadanie 4 (0.5 punktu)

Stwórz ColumnTransformer, wykonujący w pipeline'ie dla zmiennych kategorycznych (analogicznie jak w zadaniu 1 dla zmiennych numerycznych):

Dla uniknięcia dummy variable trap, oraz zmniejszenia liczby wynikowych cech, podaj drop="first" w one-hot encodingu.

Przyda się RareLabelEncoder z Feature-engine. Żeby zadziałał dla każdej cechy, podaj mu n_categories=0.

from feature_engine.encoding import RareLabelEncoder
from sklearn.preprocessing import OneHotEncoder

labels_encoder = RareLabelEncoder(n_categories=0,
                                  tol=0.01,
                                  replace_with='Rare')

one_hot_encoder = OneHotEncoder(drop='first', sparse_output=False)

categorical_pipeline = Pipeline(steps=[
    ('rare_encoder', labels_encoder),
    ('onehot', one_hot_encoder)
])

categorical_transformer = ColumnTransformer(
    transformers=[
        ('cat', categorical_pipeline, categorical_features)
    ]
)

X_train_cat = categorical_transformer.fit_transform(X_train_raw)
X_test_cat = categorical_transformer.transform(X_test_raw)

Teraz nasze cechy są w pełni przeprocesowane i gotowe do obliczeń. Trzeba je tylko połączyć - zmodyfikuj kod poniżej w zależności od nazw twoich zmiennych.

X_train = pd.concat([X_train_num, X_train_cat], axis=1)
X_test = pd.concat([X_test_num, X_test_cat], axis=1)
X_train.shape
(2045, 175)

Regresja liniowa, ridge i LASSO

Naszym baseline'owym modelem liniowym, do którego będziemy odnosić bardziej złożone modele, będzie zwykła regresja liniowa (linear regression)

Naszymi metrykami będą zarówno pierwiastek błędu średniokwadratowego (root mean squared error, RMSE), jak i średni błąd bezwzględny (mean absolute error). Oba są w tej samej jednostce, co oryginalna zmienna zależna, ale skupiają się na innych rzeczach:

W praktyce często interesuje nas kilka metryk, a mamy dużo modeli do sprawdzanie. Żeby nie kopiować kodu, warto zrobić funkcję, obliczającą i wypisującą metryki dla danego problemu. Dodatkowo tutaj musimy odwrócić logarytm, który zastosowaliśmy na początku, i takie odwrotne przekształcenia też dobrze jest zawrzeć w takiej funkcji.

from sklearn.metrics import root_mean_squared_error, mean_absolute_error


def evaluate_model(y_train: np.ndarray, y_test: np.ndarray, y_pred_train: np.ndarray, y_pred_test: np.ndarray) -> None:
    y_train = np.expm1(y_train)
    y_test = np.expm1(y_test)

    y_pred_train = np.expm1(y_pred_train)
    y_pred_test = np.expm1(y_pred_test)

    rmse_train = root_mean_squared_error(y_train, y_pred_train)
    mae_train = mean_absolute_error(y_train, y_pred_train)

    rmse_test = root_mean_squared_error(y_test, y_pred_test)
    mae_test = mean_absolute_error(y_test, y_pred_test)

    print(f"Train RMSE {rmse_train:.2f}, MAE: {mae_train:.2f}")
    print(f"Test  RMSE {rmse_test:.2f}, MAE: {mae_test:.2f}")

Jesteśmy teraz gotowi na trening naszych modeli. Na start wytrenujmy zwykłą regresję liniową.

from sklearn.linear_model import LinearRegression

reg_linear = LinearRegression()
reg_linear.fit(X_train, y_train)

y_pred_train = reg_linear.predict(X_train)
y_pred_test = reg_linear.predict(X_test)

evaluate_model(y_train, y_test, y_pred_train, y_pred_test)
Train RMSE 17643.27, MAE: 12162.87
Test  RMSE 18695.24, MAE: 13107.76

Całkiem nieźle. Dla porównania, bez usuwania skorelowanych cech i grupowania rzadkich cech, RMSE na zbiorze testowym to ponad 21 tysięcy, a model wyraźnie przeucza. Warto wykonywać odpowiedni preprocessing danych!

Mamy tu wyraźną różnicę między RMSE i MAE. Można się temu przyjrzeć i przeanalizować rozkład rezyduów (residuals), czyli błędów y − .

Zadanie 5 (0.5 punktu)

Narysuj rozkład błędów popełnianych przez model regresji liniowej. Pamiętaj o odwróceniu logarytmu eksponentą (np.expm1). Dobierz liczbę binów histogramu tak, żeby wykres był odpowiednio czytelny. Zaznacz liniami pionowymi średnią oraz medianę. Pamiętaj o opisaniu osi Y, legendzie i tytule.

Skomentuj rozkład. Jest symetryczny? Może częściej nie doceniamy wartości domu, albo ją przeceniamy, albo może są tu jakiejś pojedyncze outliery? Czy na podstawie tych rezyduów można wytłumaczyć różnicę między MAE i RMSE?

import numpy as np
import matplotlib.pyplot as plt

y_train_orig = np.expm1(y_train)
y_pred_train_orig = np.expm1(y_pred_train)
residuals_train = y_pred_train_orig - y_train_orig

y_test_orig = np.expm1(y_test)
y_pred_test_orig = np.expm1(y_pred_test)
residuals_test = y_pred_test_orig - y_test_orig

plt.figure(figsize=(18, 12))
plt.hist(residuals_train, bins=60, alpha=0.5, label='Train Residuals')
plt.hist(residuals_test, bins=60, alpha=0.5, label='Test Residuals')
plt.axvline(np.mean(residuals_train), color='r', linestyle='--', label='Mean')
plt.axvline(np.median(residuals_train), color='g', linestyle='--', label='Median')
plt.xlabel('Residuals')
plt.ylabel('Frequency')
plt.legend()
plt.title('Distribution of Linear Regression Residuals')
plt.show()

Rozkład błędów popełnianych przez model regresji liniowej wydaje się w miarę symetryczny. Istnieją pojedyncze outliery, które wpływają na zaburzenie symetrii rozkładu.

Różnica między MAE, a RMSE wynika z tego, jak każda z tych metryk traktuje błędy. W przypadku RMSE powodują one większe zaburzenie wyniku.

Mamy pewien overfitting - co prawda niewielki, ale zawsze jest coś do wyeliminowania. Dodatkowo mamy aż 175 cech - jeżeli niektóre z nich mają bardzo niską moc predykcyjną, to redukcja ich wag (lub ich usunięcie) może poprawić nawet i wynik treningowy, i testowy, bo zredukuje szum.

Warto każdy model zapisywać w zmiennej o innej nazwie, żeby później móc dokonać ich porównania. Przykładowo, reg_linear czy reg_ridge jest lepszym pomysłem niż samo reg.

Typowo siłę regularyzacji wybiera się tak, że najpierw szacujemy rząd wielkości, sprawdzając np. siatkę [1e-3, 1e-2, 1e-1, 1, 1e1, 1e2, 1e3], a później zagęszczając ją w okolicy znalezionej optymalnej wartości. Można to oczywiście robić w kilku krokach.

Zadanie 6 (1 punkt)

Wytrenuj modele z regularyzacją L2 (ridge) oraz L1 (LASSO). Dobierz siłę regularyzacji, najpierw estymując rząd wielkości, a potem zagęszczając siatkę. Dla regresji ridge wybieraj model o optymalnym MAE (argument scoring). Zmierz czas tuningu podczas drugiego kroku, po zagęszczeniu siatki.

Wypisz optymalne hiperparametry dla obu modeli, oraz ich metryki jakości.

Przydadzą się np.linspace(), RidgeCV oraz LassoCV.

from sklearn.linear_model import RidgeCV, LassoCV
from time import time

def find_best_alpha(ridge_alphas, lasso_alphas, timer = False):
    reg_ridge = RidgeCV(alphas=ridge_alphas, scoring='neg_mean_absolute_error')
    reg_lasso = LassoCV(alphas=lasso_alphas)

    models = [reg_ridge, reg_lasso]
    best_alpha = [None, None]

    for i, model in enumerate(models):
        start_time = time()
        model.fit(X_train, y_train)
        if timer:
            end_time = time()
            print(f"{model.__class__.__name__} tuning time: {round(end_time - start_time, 2)} seconds")
        print(f"{model.__class__.__name__} model alpha: {model.alpha_}")

        best_alpha[i] = model.alpha_
        y_pred_train = model.predict(X_train)
        y_pred_test = model.predict(X_test)
        evaluate_model(y_train, y_test, y_pred_train, y_pred_test)
        print()
    print()
    return best_alpha, reg_ridge, reg_lasso

reg_values = np.logspace(-3,3,7)
alpha, reg_ridge, reg_lasso = find_best_alpha(reg_values, reg_values)
new_alphas = [np.linspace(alpha[i]/10, alpha[i]*10, 1000) for i in range(2)]

_ , reg_ridge, reg_lasso = find_best_alpha(new_alphas[0], new_alphas[1], timer = True)
RidgeCV model alpha: 1.0
Train RMSE 17577.70, MAE: 12172.97
Test  RMSE 18481.38, MAE: 12973.79

LassoCV model alpha: 0.001
Train RMSE 20229.58, MAE: 13726.79
Test  RMSE 20167.56, MAE: 13954.75


RidgeCV tuning time: 5.98 seconds
RidgeCV model alpha: 1.4081081081081084
Train RMSE 17617.65, MAE: 12196.75
Test  RMSE 18468.70, MAE: 12955.05

LassoCV tuning time: 4.1 seconds
LassoCV model alpha: 0.0001990990990990991
Train RMSE 17938.56, MAE: 12366.45
Test  RMSE 18500.57, MAE: 12942.30


Dużą zaletą modeli liniowych jest możliwość analizy ważności cech oraz ich kierunku wpływu. W przypadku regresji z regularyzacją jest to szczególnie ciekawe.

Jeżeli cecha "przetrwa" z dużą wagą w ridge regression, to znaczy, że musi być naprawdę ważna, bo inaczej opłacałoby się zmniejszyć jej wagę, żeby zmniejszyć koszt z regularyzacji. Z drugiej strony, w regresji LASSO jeżeli cecha została wyeliminowana, to musiała mieć na tyle małą moc predykcyjną, że modelowi bardziej opłacało się ją usunąć i zmniejszyć koszt z regularyzacji.

Zadanie 7 (0.5 punktu)

W tym zadaniu skorzystaj z optymalnych modeli, wyznaczonych w poprzednim zadaniu.

  1. Wybierz 20 najważniejszych cech według ridge regression (w sensie wartości bezwzględnej ich wag). Przedstaw je na poziomym wykresie słupkowym (horizontal bar plot), posortowane według wartości bezwzględnej. Na wykresie mają być wartości wag z wpływem (dodatni/ujemny). Pamiętaj o opisaniu osi nazwami cech i tytule wykresu.
  2. Wypisz nazwy cech wyeliminowanych przez LASSO regression, posortowane alfabetycznie. Zakładamy, że cecha jest wyeliminowana, jeżeli dla np.isclose() jej waga jest bliska zeru.

Mogą się przydać atrybuty modelu .coef_ oraz .feature_names_in_. Wszystkie operacje można łatwo wykonać, konwertując dane na pd.Series.

T = sorted(enumerate(reg_ridge.coef_), key = lambda x: abs(x[1]), reverse=True)
x = [reg_ridge.feature_names_in_[i] for i,_ in T[:20]]
y = [j for _,j in T[:20]]
plt.bar(x, y)
plt.xticks(rotation=90)
plt.title('20 most important features by Ridge Regression')
plt.xlabel('Weight value')
plt.ylabel('Feature name')
plt.show()

lasso_coef = pd.Series(reg_lasso.coef_, index=X_train.columns)
eliminated_features = lasso_coef[np.isclose(lasso_coef, 0)].index.sort_values()
print("Features eliminated by Lasso Regression:")
print(eliminated_features)

Features eliminated by Lasso Regression:
Index(['Alley', 'BsmtCond', 'BsmtHalfBath', 'BsmtUnfSF', 'KitchenAbvGr',
       'LowQualFinSF', 'MiscVal', 'PoolArea', 'Street', 'Utilities',
       'X2ndFlrSF', 'X3SsnPorch', 'cat__BldgType_2fmCon',
       'cat__Condition1_Feedr', 'cat__Condition1_RRAe', 'cat__Condition1_RRAn',
       'cat__Condition2_Rare', 'cat__Electrical_Rare',
       'cat__Exterior1st_CemntBd', 'cat__Exterior1st_HdBoard',
       'cat__Exterior1st_MetalSd', 'cat__Exterior1st_Plywood',
       'cat__Exterior1st_VinylSd', 'cat__Exterior1st_WdShing',
       'cat__Exterior2nd_VinylSd', 'cat__Exterior2nd_Wd Shng', 'cat__Fence_No',
       'cat__Fence_Rare', 'cat__Foundation_Rare', 'cat__Foundation_Slab',
       'cat__GarageFinish_No', 'cat__GarageType_Basment',
       'cat__GarageType_BuiltIn', 'cat__GarageType_No',
       'cat__HouseStyle_1Story', 'cat__HouseStyle_2Story',
       'cat__HouseStyle_SFoyer', 'cat__LandContour_Lvl', 'cat__LotConfig_Rare',
       'cat__MSSubClass_SC120', 'cat__MSSubClass_SC190',
       'cat__MSSubClass_SC60', 'cat__MSSubClass_SC80', 'cat__MiscFeature_Rare',
       'cat__MiscFeature_Shed', 'cat__MoSold_Dec', 'cat__MoSold_Mar',
       'cat__MoSold_Oct', 'cat__Neighborhood_BrDale',
       'cat__Neighborhood_Mitchel', 'cat__Neighborhood_NAmes',
       'cat__Neighborhood_SWISU', 'cat__Neighborhood_Sawyer',
       'cat__Neighborhood_Timber', 'cat__SaleType_WD '],
      dtype='object')

Robust regression

Powyższe algorytmy nie brały pod uwagę skośnego rozkładu zmiennej zależnej, oraz realnego zapotrzebowania na medianę ceny domu. Wypróbujmy zatem algorytmy, które powinny być bardziej odporne na outliery.

Na start wypróbujemy Least Absolute Deviations (LAD) regression. Jako że jest to przypadek szczególny regresji kwantylowej, to w Scikit-learn używa się do tego klasy QuantileRegressor z podaniem odpowiedniego kwantyla.

from sklearn.linear_model import QuantileRegressor

reg_lad = QuantileRegressor(quantile=0.5)
reg_lad.fit(X_train, y_train)

y_pred_train = reg_lad.predict(X_train)
y_pred_test = reg_lad.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train, y_pred_test)
Train RMSE 80520.69, MAE: 55705.73
Test  RMSE 80920.80, MAE: 55734.18

Błąd jest wręcz tragiczny! W takich przypadkach zawsze na start trzeba się przyjrzeć predykcjom.

y_pred_test[:30]
array([11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
       11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
       11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
       11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
       11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779,
       11.99535779, 11.99535779, 11.99535779, 11.99535779, 11.99535779])

Predykcje stałe - nic dziwnego, że błąd jest duży. Scikit-learn domyślnie stosuje lekką regularyzację dla wielu modeli, ale wartość ta może być zbyt duża.

Zadanie 8 (1 punkt)

Dokonaj tuningu siły regularyzacji dla LAD regression, sprawdzając 100 wartości z zakresu od 0 do 1e-3. Wybierz model o najlepszym MAE z pomocą 5-fold CV. Zmierz czas tuningu, wypisz znalezione hiperparametry oraz sprawdź jakość finalnego modelu.

Skomentuj różnicę czasową i jakościową względem modeli ridge i LASSO. Czy w tym przypadku warto wykorzystać regresję LAD?

Przyda się GridSearchCV oraz parametr n_jobs.

from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer

alpha_values = np.linspace(0, 1e-3, 100)
reg_lad = QuantileRegressor()
param_grid = {
  "quantile" : [0.5],
  "alpha" : alpha_values}
mae_scorer = make_scorer(mean_absolute_error, greater_is_better=False)

start_time = time()
grid_search = GridSearchCV(reg_lad, param_grid, scoring=mae_scorer, cv=5, n_jobs=-1)
grid_search.fit(X_train, y_train)
end_time = time()

print(f"LAD model alpha: {grid_search.best_params_['alpha']}")
print(f"LAD tuning time: {round(end_time - start_time, 2)} seconds")

y_pred_train_lad = grid_search.predict(X_train)
y_pred_test_lad = grid_search.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train_lad, y_pred_test_lad)
LAD model alpha: 0.00047474747474747476
LAD tuning time: 642.73 seconds
Train RMSE 18170.97, MAE: 11649.58
Test  RMSE 18189.07, MAE: 12371.05

Wyniki LAD Regression są zbliżone do Ridge i LASSO, natomiast czas wykonania dużo dłuższy z powodu konieczności obliczenia układów równań programowaniem liniowym. Zatem w tym przypadku nie warto korzystać z LAD Regression.

Zamiast używać regresji LAD, opartej o programowanie liniowe, można też wykorzystać Huber regression. Jest ono mało czuła zarówno na ϵ (punkt przejścia z MSE do MAE), jak i α (siłę regularyzacji L2), więc często nie wymaga tuningu hiperparametrów.

Zadanie 9 (0.5 punktu)

Wytrenuj model Huber regression. Sprawdź jakość modelu. Skomentuj różnicę czasową i jakościową względem LAD regression (tutaj oczywiście bez tuningu). Z którego modelu skorzystałbyś w praktyce?

Domyślna liczba iteracji solwera, 100, jest zbyt mała dla naszej liczby cech, i będą z tego wynikać ostrzeżenia o braku zbieżności, więc zwiększ ją do 2000.

from sklearn.linear_model import HuberRegressor

reg_huber = HuberRegressor(max_iter=2000)
start_time = time()
reg_huber.fit(X_train, y_train)
end_time = time()

print(f"Huber tuning time: {round(end_time - start_time, 2)} seconds")

y_pred_train_huber = reg_huber.predict(X_train)
y_pred_test_huber = reg_huber.predict(X_test)
evaluate_model(y_train, y_test, y_pred_train_huber, y_pred_test_huber)
Huber tuning time: 2.34 seconds
Train RMSE 18582.51, MAE: 11701.59
Test  RMSE 17994.72, MAE: 12252.40

Huber Regression daje porównywalne wyniki (może nawet trochę lepsze), co Ridge i LASSO w stosunkowo podobnym czasie. Na plus dodatkowo jest to, że nie zwraca tak uwagi na outliery.

Use case, połączenie w system end-to-end

img.jpg

Na rynku działa firma oferująca serwis do sprzedaży nieruchomości. CEO był na prezentacji o AI, zaczął używać ChatGPT i uznał, że też musicie mieć taki potężny ML. Zostałeś zatrudniony jako (pierwszy i jedyny) data scientist, zebrałeś historyczne dane, i nawet wdrożyłeś pierwszy model regresji liniowej do estymacji wartości nieruchomości.

Pewnego dnia Product Owner przychodzi do ciebie z listą "absolutnie krytycznych" wymagań (propozycja nie do odrzucenia):

  1. Ciągle coś się sypie na produkcji, bo obecny pipeline do przetwarzania cech jest złożony z wielu kroków, które są copy-paste z twojego Jupyter Notebooka. Trzeba z tego zrobić jeden obiekt, zapisywać i wdrażać end-to-end.
  2. Nowe domy czasem mają ekstremalne wartości cech, których nie było w zbiorze treningowym, i na tych przypadkach model radzi sobie wystarczająco słabo. Może wystarczyłoby je jakoś przyciąć?
  3. Obecny model regresji liniowej działa w miarę ok, ale trzeba czegoś lepszego, co dawałoby lepsze MAE, bo większość klientów jest zainteresowana przeciętnie drogimi nieruchomościami.
  4. Sama przewidywana wartość nieruchomości to za mało, bo każdy wie, że ceny są zróżnicowane. Twoja regresja liniowa podaje tylko warunkową średnią, a zarówno sprzedający, jak i kupujący chcieliby znać sensowne górne i dolne widełki. Product Owner ze swojego doświadczenia sugeruje, że typowo takie wahania różnią się co najwyżej o 15 punktów procentowych w górę i dół od średniej ceny.

Zadanie 9 (4 punkty)

  1. Połącz wszystkie wcześniejsze transformery dla zmiennych numerycznych i kategorycznych w jeden duży ColumnTransformer. Ma do niego wejść cale X_train_raw, a wyjść gotowy output.
  2. Przycinanie wartości do zakresu znanego ze zbioru treningowego nazywa się czasem winsoryzacją (winsorization). Dodaj taki etap po skalowaniu min-max, a przed selekcją według korelacji. Przydatny będzie tutaj Feature-engine.
  3. Na podstawie wcześniejszych eksperymentów wybierz model do przewidywania własności domów. Wybór uzasadnij w komentarzu. W razie potrzeby dokonaj tuningu hiperparametrów. Zmierz jakość modelu.
  4. Wybierz i wytrenuj odpowiednie modele do regresji, żeby dostać estymację 35% i 65% ceny, w dodatku do przeciętnej ceny (estymowanej przez model z poprzedniego punktu). Dokonaj w razie potrzeby tuningu hiperparametrów. Metryką jakości może być tutaj D^2 pinball score, który ma wygodny zakres wartości [0,1] (dla dobrych modeli, im bliżej 1, tym lepiej). Podaj jakość finalnych modeli.
  5. Zapisz transformer do preprocessingu oraz estymatory do regresji do plików z pomocą Jobliba: ames_transformer.joblib, price_estimator.joblib, low_price_estimator.joblib, high_price_estimator.joblib.

Poniżej przygotowano funkcję testową, która ładuje transformer i estymatory, oraz oblicza predykcje dla trzech przykładowych domów ze zbioru testowego (taniego, przeciętnego i drogiego).

Skomentuj - czy twoim zdaniem finalne modele, podające przeciętną wartość i widełki, są subiektywnie dobrej jakości?

from feature_engine.outliers import Winsorizer
import joblib

capper = Winsorizer(capping_method='gaussian',
                    tail = "right",
                    fold=3)

num_overall_pipeline = Pipeline([("imput", median_imputer),
                                 ("scaling", min_max_scaler),
                                 ("winsorization", capper),
                                 ("selection", smart_correlated_selector)])

category_pipeline = Pipeline([("rare_labels", labels_encoder),
                              ("one_hot_encoder", one_hot_encoder)])

overall_transformer = ColumnTransformer([("numerical transform", num_overall_pipeline, numerical_features),
                                         ("category transform", category_pipeline, categorical_features)],
                                         verbose_feature_names_out=False)

X_train_overall = overall_transformer.fit_transform(X_train_raw)
X_test_overall = overall_transformer.transform(X_test_raw)

joblib.dump(overall_transformer, "ames_transformer.joblib")
['ames_transformer.joblib']
huber_true_price = HuberRegressor(max_iter=2000)
huber_true_price.fit(X_train_overall, y_train)

huber_low_price = HuberRegressor(max_iter=2000)
huber_low_price.fit(X_train_overall, np.log1p( np.expm1(y_train) * 0.35))

huber_high_price = HuberRegressor(max_iter=2000)
huber_high_price.fit(X_train_overall, np.log1p( np.expm1(y_train) * 0.65))

joblib.dump(huber_true_price, "price_estimator.joblib")
joblib.dump(huber_low_price, "low_price_estimator.joblib")
joblib.dump(huber_high_price, "high_price_estimator.joblib")
['high_price_estimator.joblib']
import joblib


def test_models() -> None:
    transformer = joblib.load("ames_transformer.joblib")
    price_reg = joblib.load("price_estimator.joblib")
    low_price_reg = joblib.load("low_price_estimator.joblib")
    high_price_reg = joblib.load("high_price_estimator.joblib")

    # select houses from test set around 25th, 50th and 75th percentile
    prices = np.expm1(y_test).astype(int)

    for price_type, quantile in [("Cheap", 0.25), ("Average", 0.5), ("Expensive", 0.75)]:
        idx = prices[prices == prices.quantile(quantile, interpolation="nearest")].index[0]
        X = X_test_raw.loc[idx].to_frame().T
        y = prices[idx]

        float_to_int = lambda x: int(np.round(np.expm1(x.item())))

        X_transformed = transformer.transform(X)
        y_pred = float_to_int(price_reg.predict(X_transformed))
        y_low_pred = float_to_int(low_price_reg.predict(X_transformed))
        y_high_pred = float_to_int(high_price_reg.predict(X_transformed))

        print(f"{price_type} house")
        print(f"True price: {y}$")
        print(f"Estimated price: {y_pred}")
        print(f"Price brackets: {y_low_pred} - {y_high_pred}")

        print("Features:")
        with pd.option_context("display.max_columns", None):
            display(X)

test_models()
Cheap house
True price: 127999$
Estimated price: 133395
Price brackets: 46691 - 86708
Features:
MSSubClass MSZoning LotFrontage LotArea Street Alley LotShape LandContour Utilities LotConfig LandSlope Neighborhood Condition1 Condition2 BldgType HouseStyle OverallQual OverallCond YearBuilt YearRemodAdd RoofStyle RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2 BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir Electrical X1stFlrSF X2ndFlrSF LowQualFinSF GrLivArea BsmtFullBath BsmtHalfBath FullBath HalfBath BedroomAbvGr KitchenAbvGr KitchenQual TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType GarageYrBlt GarageFinish GarageCars GarageArea GarageQual GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch X3SsnPorch ScreenPorch PoolArea PoolQC Fence MiscFeature MiscVal MoSold YrSold SaleType SaleCondition
195 SC50 RM 50.0 6000 2 0 4 Lvl 4 Inside 3 BrkSide Norm Norm 1Fam 1.5Fin 6 6 1940 1950 Gable CompShg MetalSd MetalSd None 0.0 3 4 CBlock 3 3 0 2 264.0 1 0.0 308.0 572.0 GasA 5 Y FuseA 848 348 0 1196 0.0 1.0 1 1 3 1 3 6 8 2 4 Detchd 1973.0 Unf 2.0 576.0 3 3 2 0 0 0 0 0 0 0 No No 0 Mar 2010 WD Normal
Average house
True price: 157999$
Estimated price: 151862
Price brackets: 53153 - 98714
Features:
MSSubClass MSZoning LotFrontage LotArea Street Alley LotShape LandContour Utilities LotConfig LandSlope Neighborhood Condition1 Condition2 BldgType HouseStyle OverallQual OverallCond YearBuilt YearRemodAdd RoofStyle RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2 BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir Electrical X1stFlrSF X2ndFlrSF LowQualFinSF GrLivArea BsmtFullBath BsmtHalfBath FullBath HalfBath BedroomAbvGr KitchenAbvGr KitchenQual TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType GarageYrBlt GarageFinish GarageCars GarageArea GarageQual GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch X3SsnPorch ScreenPorch PoolArea PoolQC Fence MiscFeature MiscVal MoSold YrSold SaleType SaleCondition
2538 SC80 RL 80.0 9600 2 0 4 Lvl 4 Inside 3 NAmes Norm Norm 1Fam SLvl 5 7 1967 1967 Gable CompShg MetalSd MetalSd BrkFace 140.0 3 3 PConc 3 3 2 5 602.0 3 402.0 137.0 1141.0 GasA 4 Y SBrkr 1141 0 0 1141 1.0 0.0 1 0 3 1 3 6 8 0 0 Attchd 1967.0 Unf 1.0 568.0 3 3 2 0 78 0 0 0 0 0 No No 0 Jul 2006 WD Normal
Expensive house
True price: 209499$
Estimated price: 190714
Price brackets: 66749 - 123940
Features:
MSSubClass MSZoning LotFrontage LotArea Street Alley LotShape LandContour Utilities LotConfig LandSlope Neighborhood Condition1 Condition2 BldgType HouseStyle OverallQual OverallCond YearBuilt YearRemodAdd RoofStyle RoofMatl Exterior1st Exterior2nd MasVnrType MasVnrArea ExterQual ExterCond Foundation BsmtQual BsmtCond BsmtExposure BsmtFinType1 BsmtFinSF1 BsmtFinType2 BsmtFinSF2 BsmtUnfSF TotalBsmtSF Heating HeatingQC CentralAir Electrical X1stFlrSF X2ndFlrSF LowQualFinSF GrLivArea BsmtFullBath BsmtHalfBath FullBath HalfBath BedroomAbvGr KitchenAbvGr KitchenQual TotRmsAbvGrd Functional Fireplaces FireplaceQu GarageType GarageYrBlt GarageFinish GarageCars GarageArea GarageQual GarageCond PavedDrive WoodDeckSF OpenPorchSF EnclosedPorch X3SsnPorch ScreenPorch PoolArea PoolQC Fence MiscFeature MiscVal MoSold YrSold SaleType SaleCondition
192 SC75 RL 0.0 7793 2 0 3 Bnk 4 Corner 3 BrkSide Norm Norm 1Fam 2.5Unf 7 7 1922 2005 Gable CompShg Wd Sdng Wd Sdng None 0.0 3 3 BrkTil 4 3 0 4 474.0 1 0.0 634.0 1108.0 GasA 3 N FuseA 1160 908 0 2068 0.0 0.0 1 1 3 1 4 8 8 1 4 Detchd 1928.0 Unf 1.0 315.0 3 3 2 0 0 60 0 0 0 0 No No 0 May 2010 WD Normal

Wybrałem Huber Regression, ponieważ nie zwraca zbyt dużej uwagi na outliery, a modele czasami niedoceniają drogich domów lub przeceniają tanie.

Zadanie dodatkowe (3 punkty)

Product Owner przychodzi do ciebie ponownie. Klienci skarżą się, że nie rozumieją, czemu model przewiduje taką, a nie inną wartość dla ich domu. Przykładowo, jeżeli to ogólny stan domu (OverallCond) czy garażu (GarageCond) obniżają cenę, to byliby gotowi je wyremontować dla lepszego zysku.

Jest to problem wyjaśnialnego AI (Explainable AI, XAI). O ile cały model w przypadku regresji liniowej jest bardzo transparentny dzięki wagom cech i liniowym predykcjom, to już zrozumienie, czemu dana predykcja była taka, a nie inna, to już problem lokalnej wyjaśnialności (local explainability), gdzie tłumaczymy pojedynczą predykcję, dla danego zestawu cech.

Najpopularniejszym algorytmem jest tutaj SHapley Additive Explanations (SHAP). Opiera się on o tzw. wartości Shapleya, mające bardzo solidne podstawy w kooperatywnej teorii gier. Opiera się na idei, że wartości cech grają w grę: jedne obniżają predykcję, a drugie podwyższają, a finalna predykcja modelu to ich konsekwencja. Każda cecha to gracz, ale jedna wnosi więcej, a inna mniej. Idea algorytmu SHAP to uznanie, że "wartościowi gracze" to wpływowe cechy, które mocno przyczyniły się do danej decyzji modelu.

Obliczanie wartości Shapleya w ogólnym przypadku ma złożoność eksponencjalną. Algorytm SHAP zaproponował efektywne ich przybliżanie dla dowolnych algorytmów ML. Dla niektórych klas modeli, przede wszystkim liniowych oraz drzewiastych, istnieją szczególnie efektywne metody, które w czasie wielomianowym obliczają dokładne wartości Shapleya. Wyjaśnienia, jak działa SHAP: link 1, link 2, link 3.

Z pomocą biblioteki shap oblicz wartości Shapleya dla przykładowych domów z ostatniego zadania i wyświetl je na tzw. force plot. Przyda ci się tutaj LinearExplainer (dokumentacja) oraz plots.force() (dokumentacja).

Porównaj ze sobą dwie metody estymacji: "interventional" oraz "correlation_dependent" (opisane w dokumentacji). Możesz je zaimplementować za pomocą maskers.Independent oraz maskers.Impute. Jako zbiór do porównania (tzw. background samples) użyj całego zbioru treningowego.

Skomentuj, czy wyjaśnienia mają twoim zdaniem sens i czy są intuicyjne. Która metoda estymacji daje subiektywnie lepsze wyniki?

import shap
from shap import maskers

transformer = joblib.load("ames_transformer.joblib")
price_reg = joblib.load("price_estimator.joblib")

background = X_train_overall

masker_independent = maskers.Independent(background)
masker_impute = maskers.Impute(background)

explainer_independent = shap.LinearExplainer(price_reg, masker_independent)
explainer_impute = shap.LinearExplainer(price_reg, masker_impute)
{"model_id":"10cd15c0ac244efe9428676a07414c3c","version_major":2,"version_minor":0}
shap.initjs()
shap.plots.force(explainer_independent.expected_value, shap_values_independent[0])
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
shap.initjs()
shap.plots.force(explainer_impute.expected_value, shap_values_impute[0])
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
shap.initjs()
shap.plots.force(explainer_independent.expected_value, shap_values_independent)
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
shap.initjs()
shap.plots.force(explainer_impute.expected_value, shap_values_impute)
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Wyjaśnienia wydają się być dość intuicyjne. Metoda Correlation Dependent wydaje się być lepsza w tym przypadku, ponieważ zakłada korelację cech i uwzględnia w obliczeniach. Metoda Interventional zakłada ortogonalność cech przez co może przekłamywać wyniki.